Numerical Simulation of the Hydrodynamical Combustion to Strange Quark Matter 
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We present results from a numerical solution to the burning of neutron matter inside a cold 
neutron star into stable u,d,s quark matter. Our method solves hydrodynamical flow equations in 
ID with neutrino emission from weak equilibrating reactions, and strange quark diffusion across the 
burning front. We also include entropy change due to heat released in forming the stable quark 
phase. Our numerical results suggest burning front laminar speeds of 0.002 — 0.04 times the speed 
of light, much faster than previous estimates derived using only a reactive-diffusive description. 
Analytic solutions to hydrodynamical jump conditions with a temperature dependent equation of 
state agree very well with our numerical findings for fluid velocities. The most important effect 
of neutrino cooling is that the conversion front stalls at lower density (below « 2 times saturation 
density). In a 2-dimensional setting, such rapid speeds and neutrino cooling may allow for a flame 
wrinkle instability to develop, possibly leading to detonation. 

PACS numbers: 97.60Jd, 26.60-c, 25.75Nq 



I. INTRODUCTION 

On grounds of asymptotic freedom in Quantum Chro- 
modynamics (QCD), hadronic matter subjected to high 
densities and/or temperatures will deconfine into a 
quark-gluon plasma. The low-density, high-temperature 
phase transition happened "in reverse" moments after 
the Big Bang, and has been fieetingly seen in ultra- 
relativistic heavy-ion collision experiments (see [l| for a 
review) . The high-density, low-temperature regime is rel- 
evant to compact stars. We assume the Witten hypoth- 
esis @: bulk strange quark matter (henceforth SQM) is 
more stable than the nuclear world we live in. The long 
lifetime of nuclei is reconciled as the improbability of ss A 
weak reactions to occur simultaneously in a nuclear vol- 
ume containing A nucleons, but SQM can still exist in the 
form of strangelets or strange quark stars, and co-exist 
with Neutron stars 0. Once SQM is nucleated inside a 
neutron star, how does it grow to form a strange quark 
star? In this paper, we numerically investigate the issue 
of combustion of pure neutron matter to u,d,s matter us- 
ing hydrodynamics, taking into account binding energy 
release and neutrino emission across the burning front - 
going beyond previous treatments of the problem 0-(3- 
This problem is interesting for two main reasons: (i) re- 
cent work [l(| shows that turbulent effects can increase 
the front velocity well beyond that expected from lami- 
nar flow analysis, entering the distributed regime which 
is a platform for subsequent detonation and (ii) conver- 
sion of a neutron star to a strange quark star has been 
investigated as an astrophysical model for gamma-ray 
bursts [TT| - F]~5| . In this work, we present an improved 
prescription of the burning front in the laminar flow ap- 
proximation, and already find speeds as high as ~ c/100, 
where c denotes speed of light. This indicates that 
unavoidable turbulent effects (such as those discussed 



in [T(| ) may well decide the fate of the conversion (defla- 
garation or detonation). Consider the situation where a 
compact star's central density has reached that of nuclear 
deconfinement, and SQM is seeded by one of many possi- 
ble alternatives Hi . Recent studies investigated the con- 
sequences of such a transition occurring during the core- 
collapse phase of a supernova [13, 13' or i ^ nucleation is 
delayed, in an older neutron star whose central density 
has increased due to spin-down [l9| . The conversion sce- 
nario we consider is non-premixed combustion (20j in a 
cold neutron star, where SQM (ash) initially grows from 
a seed by diffusion of strange quarks into neutron matter, 
viewed as a uniform (udd) mixture (fuel). The interface 
region attempts to equilibrate chemically by producing 
more strange quarks. Such a reactive-diffusive setup, as- 
suming a constant-temperature zero-thickness interface, 
was first explored in Q. Here, we consider the case for 
a macroscopically thick interface, evolved with hydro- 
dynamics, paying attention to the temperature gradient 
and neutrino emission. We find that a self-consistent nu- 
merical treatment increases the front velocity by 5-6 or- 
ders of magnitude over earlier analytic treatments [U, [2l[ . 
This large difference is mostly due to two assumptions 
made in previous analytic treatments: (i) considering the 
fluid and combustion speeds as equivalent, and (ii) lin- 
earization of the number density difference n<j — n s in 
the d + u <-> u + s reaction rate. Combustion inside a 
fluid involves flame propagation in most cases, requir- 
ing a hydrodynamical approach [9]. In addition to the 
usual conservation equations for the energy-momentum 
tensor, baryon number and electric charge, we also in- 
clude a diffusion timescale for s-quarks, neutrino emission 
and entropy evolution due to change in internal energy 
from converting to SQM. In a typical combustion, local 
temperature increase and subsequent thermal diffusion 
controls the burning rate. However, in our situation, the 
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thermal conductivity is small enough [22], |23[ that over 
the simulation time, the temperature gradient across the 
interface is unchanged. Surprisingly, this temperature 
variation becomes important through its effect on the 
pressure, not just reaction rates. We do not include dis- 
sipative terms in the hydrodynamical equations. 



II. HYDRODYNAMICS 

The 1-D hydrodynamical equations in our case are [24| : 
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and corresponding advective-diffusive terms 
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T\-§ are reaction rates for processes in Eqs.(|!|])-(13) while 
index i in the entropy source term ranges over all the par- 
ticles in the system i = {u, d, s, e~ , v}. Evolving entropy 
density s, rather than energy density, with a source term 
describing change in particle species (energy cost of "as- 
sembling" (u, d, s)-matter), allows the binding energy of 
SQM to be self-consistently taken into account. The en- 
thalpy, h, is convenient for fluids that are at relativistic 
densities. The fluid velocity, v, is expressed in units of 
the speed of light. To solve this system numerically, we 
require a constitutive equation (EoS) and the following 
reactive-diffusive inputs. 

Equation of State: In this work we use the finite- 
temperature bag model P — h — B for the EoS of SQM 
(neglecting the small electron pressure), 
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The index / in the above expressions indicates quark fla- 
vor (u,d,s). The same EoS is used for both the upstream 
(unburnt) and downstream (burnt) fluids, with the dif- 
ference being that the upstream fluid is cold u,d matter, 
since at the point of burning, the neutrons are taken to 
be already dissolved into a u,d fluid (electron pressure is 
included). We will take up the case of a more complicated 
EoS, including mixed phases, in subsequent work. 

Diffusion and Reactions: Transport of d (fuel) and s 
(ash) quarks through the interface driven by concentra- 
tion gradients results in colliding flows of different flavors. 
The diffusion coefficient relevant for burning into SQM 
is 
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Equilibrium in SQM is established by beta-decay and 
electron capture reactions, 
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We use rates given by |25| , which also have equilibrium- 
seeking terms for the leptonic processes, 
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where Afj, = (fi d — fj, s ), pp is the quark's Fermi mo- 
mentum, Gp Fermi's constant and 9c the Cabbibo angle. 
For process (|13p to proceed in a given region, a mini- 
mum number of s quarks must be present. While this 
number should depend on factors such as the strangelet 
mass and surface tension, here we simply make sure to 
avoid unphysical effects, such as superluminous diffusion 
speeds [24] , by imposing a smooth cut-off on the s quark 
Fermi momentum {pf s = \/ H 2 — fn 2 ^ 0.1 MeV) for reac- 
tion (|13|) to proceed (this is analogous to the activation 
temperature, in Arrhenius-type reactions, typically used 
in modeling heat-diffusion driven combustion). 

Neutrino emission: Neutrinos are emitted copiously 
from the location of the interface, where the leptonic 
weak reaction rates from chemical equilibration are high- 
est. At these temperatures (tens of MeV) and densities 
(p ~ 10 15 g/cc), neutrino mean free paths A are on the or- 
der of 100 cm [261 ] . Accurate neutrino transport requires 
solving the Boltzmann equation, which in our setup intro- 
duces additional stiffness in the flow equations. A simpler 
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estimate capturing the essential physics in 1-D is to in- 
troduce an exponential cut-off on the neutrino emissivity 
as follows 

e = {e qf) + e qs )xe-^- x ^ x . (18) 

where e q p and e qs denote the non-equilibrium neutrino 
emission rate for reactions © and (|TTj) respectively [25[ , 
x is the position of the emitting region and x\ is the po- 
sition of the front of the burning interface. Effectively, 
for a given emitting region at x, if the the distance to 
the interface x\ — x is more than the mean free path A, 
then produced neutrinos are trapped, otherwise they es- 
cape. Since A ~ 100cm, neutrinos produced near or in the 
interface and directed outwards essentially free stream. 
The matter ahead of the burning interface is cool, while 
SQM behind it is hot and produces many neutrinos. The 
small mean free path then implies that neutrino cool- 
ing does not significantly alter the temperature of equili- 
brated SQM on the timescale of the simulation, but has 
an important effect on the diffusion of strange quarks 
across the interface, and hence the speed of the burning 
front (Fig. El). 

III. NUMERICAL SIMULATIONS & RESULTS 

The variables in the equations of hydrodynamical com- 
bustion (Eq. are solved for numerically using a fourth- 
order Runge-Kutta scheme. Spatially, a third-order up- 
winded advection, flux-limited, finite- volume approach is 
used [27|] . The diffusion and pressure gradient terms are 
second-order, not upwinded, and treated separately from 
the advection terms (ie. not flux- limited). A large pres- 
sure wave is created from the initial state, even though 
it initially satisfies pressure equilibrium. This is typical 
for combustion problems, where the (unburnt) fluid in 
front of the interface is set in motion (|24|, pg. 487). 
However, this wave is transient and quickly flows past 
the burning region, increasing the speed of the inter- 
face without impacting its long-term evolution (for ani- 
mations, visit http:/ /www.capca.ucalgary.ca) An accept- 
able grid spacing is Ax = 0.05, resulting in a limited 
timestep of At/ Ax < 0.3 from the advection terms, and 
At/ {Ax) z < 1/D from the diffusion terms. Leaving 
more detailed description of numerical aspects to a sub- 
sequent article, we discuss here our main physical results. 

(i) Effects of hydrodynamics: In Fig. Q] the interface 
speed, with and without the effects of hydrodynamics, 
is plotted for various initial conversion densities. In the 
former case, typical speeds for the burning interface were 
found to be between 0.002c and 0.04c for initial baryon 
densities ranging from 1.7no to 5.3no> where uq is nu- 
clear saturation density. These burning speeds are much 
higher than previous estimates [3, [2l| . The reasons come 
from the T and /i s variations across the finite-width in- 
terface. Just after contamination, at small values of fi s , 
the reactions producing s quarks are dominated by the 



Afi 3 factor in Eqn. (lT7"l) . Further behind the interface, 
s quark production becomes increasingly dependent on 
the temperature term. This increases the reaction rate 
as expected [28| . resulting in a faster speed of the burn- 
ing interface. Including hydrodynamics in the reactive- 
diffusive simulations creates different fluid velocities on 
either side of the interface, which ends up effectively 
opposing the interface's progression (discussed below). 
Typical widths of the interface (see Fig. [2]) were found to 
be ~ 1 cm when hydrodynamics is included and ~ 10 cm 
for a purely reaction-diffusion system. 
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FIG. 1. Steady-state burning interface speeds Wburn for sim- 
ulations with various initial densities (quark chemical po- 
tential) /iiNiT- The three hydrodynamic cases (HYDRO) 
are with no neutrino cooling (dashed line), neutrino cooling 
from Eg 1181 (solid line), and enhanced neutrino cooling (dash- 
dotted line). The dotted line indicates simulations without 
hydrodynamics (ie. fluid velocities zero everywhere). Ubm-n 
increases with larger densities, since more fuel is present, and 
decreases with larger cooling rates. As cooling becomes more 
effective, the hydrodynamic jump conditions (Eq, IAl[l are sat- 
isfied by increasingly opposing the advance of the interface, 
which consequently stalls at progressively higher densities. 

(ii) Effects of z^-cooling: Neutrino emission (delep- 
tonization) causes a decrease in pressure for the burnt 
fluid. The resulting pressure gradient forces fluid ve- 
locities to become increasingly negative (in the refer- 
ence frame comoving with the burning interface), caus- 
ing advection to oppose the progression of the burning 
interface [29j . Since cooling rates may have uncertain- 
ties, we parameterize the efficacy of neutrino cooling by 
C = T — Tcooiod, where T coo i c d and T are the down- 
stream (burnt) temperatures with and without cooling. 
As shown by the two temperature profiles in Fig. [5J Even 
a modest drop in temperature can decrease the pressure 
enough to enter an advection dominated regime, where 
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the upstream fluid velocity (vi) advects the interface 
backwards faster than it can progress due to reactions 
and diffusion (|ui| > wrd)- In such a case the interface 
halts, as seen in Fig. [3J as soon as urd + i>i < 0. This 
is because as the interface stops diffusing into the fuel, 
the reactions are no longer proceeding. Since neutrino 
production drops as a consequence, energy is no longer 
being removed from the burning region and the system 
reaches a situation where diffusion and advection are in 
balance. 



> 

CD 




0.0' 



400. MeV 
350 MeV 



300 MeV 



0.10 1.00 
C [MeV] 



0.00 



[cm] 



FIG. 3. Velocity of the burning interface, purely from the 
reaction-diffusion process (vrb\ ie. no hydrodynamics) plus 
the upstream fluid velocity (vi), versus cooling. The upstream 
velocity is calculated analytically from the jump conditions 
(cf. Apx. [XJ, and neutrino cooling is represented by the 
difference between non-cooled and cooled downstream tem- 
peratures C — T — T coo i c( j. Values shown are the initial den- 
sities. The interface halts after a critical amount of energy is 
removed by cooling. 



FIG. 2. A snapshot during the simulation of the temperature 
(T), reaction rate (R), and neutrino emissivity, e„, throughout 
the burning interface. The temperature is shown with (solid 
line) and without (dashed line) neutrino cooling effects, where 
the difference between the two is the variable C = T — T coo i e d 
that serves as the measure of cooling. 



IV. DISCUSSION & CONCLUSIONS 

We have performed I-D numerical simulations of the 
burning of neutron matter to strange quark matter 
(SQM) with consistent treatment of reactions, diffusion, 
and hydrodynamics. By modeling a region of SQM sur- 
rounded by u,d matter, interpreted as neutron matter 
above nuclear densities, we find typical speeds of the 
burning process to be between 0.002c and 0.04c and in- 
terface widths of ~ 1 cm. These speeds are noticeably 
hi gher than estimates found in previous works, for eg., 
[til |22| . In this work we have addressed the importance of 
evolving temperature self-consistently from the binding 
energy release during conversion to SQM by incorporat- 
ing this with the entropy evolution equation. We have 
also shown how neutrino cooling can halt the burning 
interface by decreasing pressure support against advec- 
tive forces. Hence, the importance of neutrinos cannot 



be overstated and must be addressed more thoroughly 
in future work. An equally important focus for future 
work is a two-dimensional treatment. While cooling can 
only halt the interface in one dimension, in two or more 
dimensions we propose that a new type of instability 
would develop, caused by regions along the burning in- 
terface halting due to cooling, at which point unburnt 
material starts to flow backwards onto the interface (as 
inferred from the jump conditions, Apx. |XJ), whereas re- 
gions not halted by cooling will have unburnt material 
flowing away from the interface. The result is a wrinkled 
interface, with shearing between the unburnt fluids of 
halted and non-halted regions. A wrinkled interface in- 
creases the diffusion rate, and causes an overall increase 
in the burning interface's speed (wbum)- However, the 
wrinkling is also subject to stabilization by diffusion - 
in the dimension along the interface concave regions are 
accelerated while convex regions are decelerated. The 
interplay between stabilization and the wrinkling insta- 
bility can result in three scenarios: either (i) stabilization 
is too strong causing Wbum to remain small and the entire 
interface halts, or (ii) stabilization is moderate and the 
interface progresses outwards as a combustion, or (iii) 
stabilization is weak and Ubum increases without bound, 
presumably resulting in a detonation. Validating these 
options would require high-resolution multi-dimensional 
simulations, which we leave for future work. 
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Appendix A: Verification of Numerical Solutions 

We performed numerical tests separately for the diffu- 
sive, reactive, and hydrodynamic parts of the code. With 
only diffusion, we run the usual tests of a diffusing ini- 
tially gaussian profile and find a relative error in the gaus- 
sian width that is smaller than the resolution at all times. 
The reactive part of the code confirmed analytically esti- 
mated timescales to achieve weak equilibrium [J, |25l |28| . 
For hydrodynamics, we solved jump conditions in the 
frame of the burning interface, including the tempera- 
ture increase due to the release of binding energy. From 
Eq.©, 
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The subscripts 1 and 2 indicate upstream (unburnt) and 
downstream (burnt) fluids respectively. Pressure is given 
in terms of enthalpy ((Eq. [5]). The above expressions 
are solved analytically, and upstream and downstream 
velocities agree to better than 2% with those found nu- 
merically. We do not include the jump condition from 
the entropy equation, since the reaction term introduces 
a non-linear component, so the temperature increase due 
to release of binding energy is not found analytically. In- 
stead, we used computed values from simulations without 
hydrodynamics. 



